

% streamlines some of the functionality of widefieldProcessingScripts.m
% 
% optional argument backgroundCorrectionType can be 'avg' or 'frame'; 'avg'
% will just subtract the average backgroudn whereas 'frame' will subtract
% the frame-by-frame background if that has been acquired
%
%
% updates as of 8/23/2023:
% removed dropped frames and hot pixel correction since that's not an issue
% with the Dalsa camera
%
% update as of 6/18/24:
% adding in an option for movtype to be 1 channel IOS (called 'ios1ch') as
% well as optogenetic stimulation
%
% update as of 6/27/2024:
% adding an option of 'whisker' as a stim type
%
% update as of 7/2/2024:
% adding normalizeSingleTrials as an input. this is a logical that if set
% to true will run the script as normal, normalizing each stim to the
% baseline but if set to false it will skip that step
function [r,mov,sortedCh1,sortedCh2,sortedCh1Z,sortedCh2Z,ontarget_gcamp,channelNormalization,baselinePixelValues] = widefieldProcessingScripts_v2_2(movname,metaname,stimtype,movtype,wavelengths,driftCorrectionTransform,dosave,mov,windowmask,backgroundCorrectionType,background,binning,gcampcorrection,normalizeSingleTrials,bgmetaname)


if nargin < 8 || isempty(mov)
    moviepreloaded = false;
else
    moviepreloaded = true;
end
if nargin < 13 || isempty(gcampcorrection)
    gcampcorrection = false;
end
if nargin < 12 || isempty(binning)
    binning = 1;
end
if nargin < 14 || isempty(normalizeSingleTrials)
    normalizeSingleTrials = true;
end
if nargin < 15 || isempty(bgmetaname)
    bgmetaname = [];
end
if nargin < 1 || isempty(movname)
    homedir = cd;
    [fn,fd] = uigetfile('*.tif','select one .tif file in the sequence');
    cd(fd)
    movname = fn;
    im = loadTiffStack(fn,1:6);
    cd(homedir)
elseif ischar(movname) 
    temp = dir([movname,'*']);
    if moviepreloaded
        im = mov(:,:,1:6);
    else
        im = loadTiffStack(temp(1).name,1:6);
    end
end

if nargin <  2 || isempty(metaname)
    metaname = uigetfile('*.mat','select metadata file');
end

if nargin < 11 || isempty(background)
    temp = dir('*ackground*.tif');
    background = temp(1).name;
end
if nargin < 10 || isempty(backgroundCorrectionType)
    backgroundCorrectionType = 'frame';
end

if nargin < 4 || isempty(movtype)
    movtype = 'gcamp';      %'gcamp' or 'ios'
end
if nargin < 5 || isempty(wavelengths)
    switch movtype
        case 'gcamp'
            wavelengths = [];       %blank for 'gcamp'; [636,530] for 'ios'
        case 'ios'
            wavelengths = [636,530];
        case 'ios1ch'
            wavelengths = [530];
    end
end
if nargin < 6 || isempty(driftCorrectionTransform)
    driftCorrectionTransform = [];      %empty if no previously saved drift correction. otherwise the filename of the drift correction file. This file is saved as an output of parsePeriodicStimMovies_v4_1
end
imagefiltering = [true,true];
donormalize = true;
if nargin < 3 || isempty(stimtype)
    stimtype = 'stationary';
end
if nargin < 7 || isempty(dosave)
    dosave = true;
end

%% select the windowmask:
if nargin < 9 || isempty(windowmask)
    figure
    im = median(im,3);
    imagesc(im);
    roi = drawpolygon;
    windowmask = createMask(roi);
end
%% load in the tif stack:
if nargin < 8 || isempty(mov)
    rootname = movname(1:(length(movname)-6));
    
    if contains(movname,'.tif')
        mov = loadTiffStack2(movname);
    else
        mov = loadTiffStack2([movname,'.tif']);
    end
end


%% run processing
if strcmp(stimtype,'optogenetic')
    [avMovie_ch1, avMovie_ch2, vascReference,avHbO, avHbR,uniqueIdentifier,~,~,sortedCh1,sortedCh2,stdMovie_ch1,stdMovie_ch2,sortedCh1Z,sortedCh2Z,channelNormalization,baselinePixelValues] = ...
        parsePeriodicStimMoviesOPTO_v4_4(mov,background,metaname,movtype,imagefiltering,gcampcorrection,driftCorrectionTransform,wavelengths,windowmask,donormalize,stimtype,backgroundCorrectionType,bgmetaname);
    load(metaname,'totalFrameTime','stimONframe','stimOFFframe','laserpower')
    fps = 1/(totalFrameTime/1000);
    preStimTime = stimONframe/fps;
    stimlength = (stimOFFframe-stimONframe)/fps;
elseif strcmp(stimtype,'whisker')
    [avMovie_ch1, avMovie_ch2, vascReference,avHbO, avHbR,uniqueIdentifier,~,~,sortedCh1,sortedCh2,stdMovie_ch1,stdMovie_ch2,sortedCh1Z,sortedCh2Z,channelNormalization,baselinePixelValues] = ...
        parsePeriodicStimMoviesWHISKER_v4_4(mov,background,metaname,movtype,imagefiltering,gcampcorrection,driftCorrectionTransform,wavelengths,windowmask,donormalize,stimtype,backgroundCorrectionType,normalizeSingleTrials,bgmetaname);
    load(metaname,'totalFrameTime','stimONframe','stimOFFframe','whiskerpower')
    fps = 1/(totalFrameTime/1000);
    preStimTime = stimONframe/fps;
    stimlength = (stimOFFframe-stimONframe)/fps;
elseif strcmp(stimtype,'moving')
    [avMovie_ch1, avMovie_ch2, vascReference,avHbO, avHbR,uniqueIdentifier,~,~,sortedCh1,sortedCh2,stdMovie_ch1,stdMovie_ch2,sortedCh1Z,sortedCh2Z,channelNormalization,baselinePixelValues] = ...
        parsePeriodicStimMovies_MOVING_v1_0(mov,background,metaname,movtype,imagefiltering,gcampcorrection,driftCorrectionTransform,wavelengths,windowmask,donormalize,stimtype,backgroundCorrectionType,normalizeSingleTrials,bgmetaname);
    load(metaname,'totalFrameTime','stimparams','preStimTime','framesPerStim','fps','rectPositionsDrift','movindicator')
    rectPositions = rectPositionsDrift;
    stimlength = stimparams{1}.totalStimLength;
else
    [avMovie_ch1, avMovie_ch2, vascReference,avHbO, avHbR,uniqueIdentifier,~,~,sortedCh1,sortedCh2,stdMovie_ch1,stdMovie_ch2,sortedCh1Z,sortedCh2Z,channelNormalization,baselinePixelValues] = ...
        parsePeriodicStimMovies_v4_4(mov,background,metaname,movtype,imagefiltering,gcampcorrection,driftCorrectionTransform,wavelengths,windowmask,donormalize,stimtype,backgroundCorrectionType,normalizeSingleTrials,bgmetaname);
    load(metaname,'totalFrameTime','stimparams','preStimTime','framesPerStim','fps','rectPositions','movindicator')
    stimlength = stimparams{1}.totalStimLength;
end
if strcmp(movtype,'gcamp')
    gcampstims = zeros(size(avMovie_ch1{1},1),size(avMovie_ch1{1},2),length(avMovie_ch1));
    iosstims = gcampstims;
    gcampstimStd = gcampstims;
    iosstimStd = gcampstims;
elseif strcmp(movtype,'ios')
    ios1stim = zeros(size(avMovie_ch1{1},1),size(avMovie_ch1{1},2),length(avMovie_ch1));
    ios2stim = zeros(size(avMovie_ch2{1},1),size(avMovie_ch2{1},2),length(avMovie_ch2));
    avHbOstims = ios1stim;
    avHbRstims = ios1stim;
    totalFrameTime = (1/fps)*2*1000;
elseif strcmp(movtype,'ios1ch')
    iosstims = zeros(size(avMovie_ch2{1},1),size(avMovie_ch2{1},2),length(avMovie_ch2));
    iosstimStd = iosstims;
end

sf1 = (round(preStimTime/(totalFrameTime/1000))+1):round((preStimTime+1)/(totalFrameTime/1000));
sf2 = (round((preStimTime+1)/(totalFrameTime/1000))+1):round((preStimTime+stimlength+.5)/(totalFrameTime/1000));
for n = 1:length(avMovie_ch1)
    if strcmp(movtype,'gcamp')
        gcampstims(:,:,n) = nanmean(avMovie_ch1{n}(:,:,sf1),3);
        iosstims(:,:,n) = nanmean(avMovie_ch2{n}(:,:,sf2),3);
        gcampstimStd(:,:,n) = nanmean(stdMovie_ch1{n}(:,:,sf1),3);
        iosstimStd(:,:,n) = nanmean(stdMovie_ch2{n}(:,:,sf2),3);
    elseif strcmp(movtype,'ios')
        ios1stim(:,:,n) = nanmean(avMovie_ch1{n}(:,:,sf2),3);
        ios2stim(:,:,n) = nanmean(avMovie_ch2{n}(:,:,sf2),3);
        avHbOstims(:,:,n) = nanmean(avHbO{n}(:,:,sf2),3);
        avHbRstims(:,:,n) = nanmean(avHbR{n}(:,:,sf2),3);
    elseif strcmp(movtype,'ios1ch')
        iosstims(:,:,n) = nanmean(avMovie_ch2{n}(:,:,sf2),3);
        iosstimStd(:,:,n) = nanmean(stdMovie_ch2{n}(:,:,sf2),3);
        avMovie_ch1 = [];stdMovie_ch1 = [];
    end
end
writetiff(uint16(vascReference),'lowResVasculatureReference01.tif');

r.avMovie_ch1 = avMovie_ch1;
r.avMovie_ch2 = avMovie_ch2;
r.stdMovie_ch1 = stdMovie_ch1;
r.stdMovie_ch2 = stdMovie_ch2;
r.vascReference = vascReference;
r.avHbO = avHbO;
r.avHbR = avHbR;
r.uniqueIdentifier = uniqueIdentifier;
r.windowmask = windowmask;

if dosave
    saveroot = movname;
    saveroot = strrep(movname,'.tif','');

    if strcmp(movtype,'gcamp')
        temp = [saveroot,'gcampAvStimsUncorrected_normalized.tif'];
        writetiff(uint16(gcampstims*1000),temp);
        temp = [saveroot,'iosAvStims_normalized.tif'];
        writetiff(uint16(iosstims*1000),temp);
        temp = [saveroot,'averageStims_uncorrected_normalized.mat'];
        save(temp,'r','sf1','sf2','rectPositions','stimparams','movindicator','uniqueIdentifier','channelNormalization','baselinePixelValues')
        temp = [saveroot,'iosStdStims.tif'];
        writetiff(uint16(iosstimStd*1000),temp);
        temp = [saveroot,'gcampStdStims.tif'];
        writetiff(uint16(gcampstimStd*1000),temp);
    elseif strcmp(movtype,'ios')
        temp = [saveroot,'averageStims_normalized.mat'];
        save(temp,'avMovie_ch1','avMovie_ch2','windowmask','sf1','sf2','avHbR','avHbO','uniqueIdentifier','vascReference','rectPositions','stimparams','movindicator','channelNormalization','baselinePixelValues')
        temp = [saveroot,'HbOavg.tif'];
        writetiff(int16(avHbOstims*10^9),temp)
        temp = [saveroot,'HbRavg.tif'];
        writetiff(int16(avHbRstims*10^9),temp)
    elseif strcmp(movtype,'ios1ch') & strcmp(stimtype,'optogenetic')
        temp = [saveroot,'averageStims_normalized.mat'];
        ontarget_gcamp = [];
        save(temp,'avMovie_ch2','windowmask','sf2','avHbR','avHbO','uniqueIdentifier','vascReference','laserpower','channelNormalization','baselinePixelValues')
    end
end

%% plot the gcamp values in the activated areas over the whole timeseries
% generate masks from the average stim images:
if ~strcmp(movtype,'ios1ch')
    thresholdMin = 3;
    masks = false(size(r.avMovie_ch1{1},1),size(r.avMovie_ch1{1},2),length(r.avMovie_ch1));
    for n = 1:size(masks,3)
        im = mean(r.avMovie_ch1{n}(:,:,sf1),3,'omitnan');
        im = im - prctile(im(:),thresholdMin);
        im(im<0) = 0;
        im = im/max(im(:));
        t = multithresh(im,3);
        masks(:,:,n) = im>t(3);
    end

    % get average gcamp value in the on taret area for each stim over each
    % timepoint
    ontarget_gcamp = nan(length(sortedCh1)*size(sortedCh1{1},4),size(sortedCh1{1},3));
    count = 0;
    for n = 1:size(masks,3)
        for n2 = 1:size(sortedCh1{n},4)
            count = count + 1;
            for n3 = 1:size(sortedCh1{n},3)
                cur = sortedCh1{n}(:,:,n3,n2);
                ontarget_gcamp(count,n3) = mean(cur(masks(:,:,n)));
            end
        end
    end
    figure
    temp = ontarget_gcamp(:,2:end);
    cl = [min(temp(:)),max(temp(:))];
    imagesc(ontarget_gcamp,cl);colormap 'jet';
    
    avggcamptraj = zeros(length(r.avMovie_ch1),size(sortedCh1{1},3));
    for n = 1:length(r.avMovie_ch1)
        idx = ((n-1)*size(sortedCh1{n},4)+1):((n)*size(sortedCh1{n},4));
        avggcamptraj(n,:) = mean(ontarget_gcamp(idx,:),1);
    end
    figure
    try
        t = (0:(size(avggcamptraj,2)-1))*totalFrameTime/1000;
        plot(t,avggcamptraj','linewidth',2)
        hold on
        yl = [.9,1.1];
        ylim(yl);
        h = fill([preStimTime,stimparams{1}.totalStimLength+preStimTime,stimparams{1}.totalStimLength+preStimTime,preStimTime],...
            [yl(1),yl(1),yl(2),yl(2)],[.5,.5,.5]);
        h.LineStyle = 'none';h.FaceAlpha = .3;
    catch
        disp('error in some plotting')
    end
end

